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We study off-equilibrium behaviors at first-order transitions (FOTs) driven by a time dependence 
of the temperature across the transition point T c , such as the linear behavior T(t)/T c = 1 ± t/t s 
where t 3 is a time scale. In particular, we investigate the possibility of nontrivial off-equilibrium 
scaling behaviors in the regime of slow changes, corresponding to large t s , analogous to those arising 
at continuous transitions, which lead to the so-called Kibble-Zurek mechanism. 

We consider the 2D Potts models which provide an ideal theoretical laboratory to investigate issues 
related to FOTs driven by thermal fluctuations. We put forward general ansatzes for off-equilibrium 
scaling behaviors around the time t = 0 corresponding to T c . Then we present numerical results 
for the q — 10 and q = 20 Potts models. We show that phenomena analogous to the Kibble- 
Zurek off-equilibrium scaling emerge also at FOTs with relaxational dynamics, when appropriate 
boundary conditions are considered, such as mixed boundary conditions favoring different phases at 
the opposite sides of the system, which enforce an interface in the system. 

PACS numbers: 05.70.Fh,64.70.qj,64.60.an 


I. INTRODUCTION 

Slow time variations of system parameters across con¬ 
tinuous transitions inevitably lead to off-equilibrium be¬ 
haviors, giving rise to the so-called Kibble-Zurek (KZ) 
mechanism P,[2] . A typical example is provided by a lin¬ 
ear change of the temperature across its transition value 
T c , such as T(t)/T c = 1 — t/t s starting from t = t t < 0 
to t = tf > 0, where t s controls the speed of the temper¬ 
ature variation. The emergence of off-equilibrium phe¬ 
nomena is essentially related to the fact that continu¬ 
ous transitions develop correlations with diverging length 
scale, which cannot adapt themselves to the variations of 
the temperature, even in the regime of slow variations. 
However, in the limit of large t s , the system develops 
an off-equilibrium scaling behavior involving t s , which is 
controlled by the same critical exponents of the system 
at equilibrium 0-i- This issue has been also extended 
to quantum transitions, obtaining analogous behaviors 
when quasi-adiabatic changes of an external parameter 
go through continuous quantum transitions [5)-j7|. Sev¬ 
eral experiments have investigated these off-equilibrium 
phenomena, in particular checking the predictions for 
the abundance of residual defects arising from the off- 
equilibrium conditions across T Cl as predicted by the KZ 
mechanism, see, e.g., Refs. IsWfSOll . 

In this paper we investigate whether analogous phe¬ 
nomena arise in systems undergoing first-order transi¬ 
tions (FOTs), characterized by a discontinuity of the en¬ 
ergy density in the thermodynamic limit. Unlike contin¬ 
uous transitions, the length scale of the correlations in 
the thermodynamic limit (within each phase) remains fi¬ 
nite when approaching T c . Nevertheless, we show that 
off-equilibrium scaling behaviors, analogous to the KZ 
scaling at continuous transitions, may also arise at FOTs 
when slowly varying the temperature across T c , in sys¬ 
tems with appropriate boundary conditions favoring the 


presence of an interface. 

Two-dimensional (2D) g-state Potts models provide an 
ideal theoretical laboratory to investigate issues related 
to FOTs driven by the temperature (when q > 4). We 
consider the off-equilibrium behavior arising from a re¬ 
laxational dynamics with a time-dependent temperature 
T crossing T c , such as T(t)/T c « 1 ±t/t s . We show that, 
when slowly crossing the FOT, i.e., for large time scales 
t s , the off-equilibrium behavior turns out to be depen¬ 
dent on the geometry and boundary conditions. This is 
essentially related to the dependence of the equilibrium 
relaxational dynamics at FOTs on the boundary condi¬ 
tions. For symmetric boundary conditions, such as peri¬ 
odic boundary conditions (PBC), systems of size L are 
characterized by an exponentially slow dynamics due to 
an exponentially large tunneling time r ~ e aL between 
the coexisting phases. On the other hand, power-law 
behaviors characterize the slow dynamics when mixed 
boundary conditions (MBC) are considered, i.e., when 
the boundary conditions at two opposite sides of the sys¬ 
tem are related to the different high-T and low-T phases, 
effectively generating an interface. We argue that the 
MBC settings lead to a power-law off-equilibrium scaling 
behavior involving the time scale t s of the slow temper¬ 
ature variation across the FOT point. This is confirmed 
by a numerical analysis of Monte Carlo (MC) simula¬ 
tions following a protocol similar to those considered for 
the KZ problem at continuous transitions. 

The paper is organized as follows. In Sec. [Tl]we present 
the 2D Potts model in which we develop and check the 
off-equilibrium scaling theory at FOTs. There we also de¬ 
fine the protocol we consider for the time variation of the 
temperature across T c , which leads to the off-equilibrium 
behavior. In Sec. m we develop an off-equilibrium scal¬ 
ing theory at FOTs, stressing the crucial dependence on 
the geometry and boundary conditions of the system un¬ 
dergoing the FOT. We essentially report results for the 
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2D Potts model, but the main features can be straight¬ 
forwardly generalized to other systems. In Sec . 11VI we re¬ 
port a numerical analysis which provides support to the 
off-equilibrium scaling behavior put forward in Sec. IIIII 
Finally, we draw our conclusions in Sec. |V] The appen¬ 
dices contain some details of our numerical study. 


II. THE MODEL AND THE OFF-EQUILIBRIUM 
PROTOCOL 


A. The 2D Potts model 


2D g-state Potts models provide a useful theoretical 
laboratory where to study issues related to FOTs driven 
by thermal fluctuations. They are defined by the parti¬ 
tion function 

Z=J2?- H/T , H = -'£,8(8 x ,Sy), (1) 

{«*} ( x y) 

where the sum in the Hamiltonian is meant over the 
nearest-neighbor sites of a square lattice, s x are integer 
variables 1 < s x < g, S(a,b) = 1 if a = b and zero 
otherwise. The Potts model undergoes a phase transi¬ 
tion [Ud^l at 


Pc = T~ l = ln(l + V?), (2) 

which is continuous for q < 4 and first order for q > 
4. For q = 2 the Potts model becomes equivalent to 
the Ising model. FOTs for q > 4 become stronger and 
stronger with increasing g, indeed the latent heat grows 
with increasing q. 

We consider 2D Potts models with two different ge¬ 
ometries: square L x L lattices and anisotropic L± x L\\ 
slab-like lattices with L\\ L±. We also consider differ¬ 

ent boundary conditions: periodic boundary conditions 
(PBC) and mixed boundary conditions (MBC) where op¬ 
posite boundary sides are related to the different high-T 
and low-T phases. 

More precisely, in the MBC case we consider 
anisotropic Lj_ x L\\ lattices with L± = 2L + 1, so that 
— L < x\ < L and 1 < Xi < Ly. We take open boundary 
conditions along the x± = L line, which corresponds to 
having T = oo bonds between the line sl, X2 and a further 
fictitious line sl+i,x 2 - At the opposite side, in order to 
have boundary conditions corresponding to the ordered 
T = 0 phase, we add a fictitious x± = —L — 1 line where 
we fix S-l- i, X2 = 1, and add the corresponding bond 
terms to the Hamiltonian. The boundary conditions are 
chosen periodic along the X 2 direction of size L\\ . Note 
that MBC breaks explicitly the g-state permutation sym¬ 
metry of the Potts model. 

In our study we consider observables related to the 


magnetization and energy density, i.e., 
1 ^ qS(s x , 1) - 1 


m 




V 

1 


q-1 


e Y ^ ' ^(^11,3:2 1 S Xl,X 2 + l)j 


( 3 ) 

( 4 ) 


where V is the number of sites of the lattice. Their 
equilibrium values at the FOT point are known for any 
q > 4 [32]. In particular, approaching the transition point 
after the thermodynamic L —> 00 limit, we have 


e~ = e(T“) = 0.910342..., e+ = e(T +) = 0.313265..., 
m c = m(T~) = 0.9411759..., for q = 20, (5) 

and 

e“ = 0.832126..., e+= 0.428553..., 

m c = 0.857106..., for q = 10. (6) 


We also define the related renormalized quantities 


-(*) = 


i(t) 


m c 


-(t) = — 






( 7 ) 


so that, at equilibrium and in the thermodynamic limit, 
m r = e r = 0 for T T+ and m r = e r = 1 for T T~ . 

The correlation length related to the exponential decay 
of the two-point function in the limit T —> Tp [after 
the thermodynamic limit) is exactly known [331.1341). It 
decreases with increasing g, e.g., £ + = 2.6955... for q = 20 
and £ + = 10.5595... for q = 10. Numerical results [35l[36| 
support the hypothesis that the correlation length £~ for 
T T~ equals £+. 


B. Off-equilibrium protocol across the transition 

The protocol that we consider for the off-equilibrium 
simulations across T c is similar to that leading to the KZ 
mechanism a a at continuous transitions. We vary the 
temperature across the transition point and study the re¬ 
sulting off-equilibrium behavior in the limit of slow time 
variations. More precisely, we vary the inverse tempera¬ 
ture (3 = 1/T so that 

6(t) = fi(t)/p c - 1 = ±t/t s (8) 

where t € [t t < 0,t/ >0] is a time variable varying from 
a negative to a positive final value. The value t = 0 cor¬ 
responds to 6(t) = 0, i.e., T(t) = T c . The parameter t B 
provides the time scale of the temperature variation. We 
start our simulations from equilibrium configurations at 
P = P c [l + <5(tj)]. Then we make the system evolve un¬ 
der a purely relaxational dynamics (also known as model 
A of critical dynamics [37|), which can be realized by 
standard Metropolis or heat bath updatings in MC sim¬ 
ulations, see App. 0 The ± sign in Eq. © corresponds 
to crossing T c starting from the high-T phase (sign +) 
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or from the low-T phase (sign —). The unit time for t 
corresponds to a complete sweep of the whole lattice by 
heat bath or Metropolis upgradings. The temperature is 
changed according to Eq. © every sweep, incrementing 
t by one. We stop the off-equilibrium relaxational dy¬ 
namics when t = tf. We repeat this procedure several 
times averaging the observables at fixed time f, thus the 
average is performed on the equilibrium Gibbs ensemble 
of the initial inverse temperature /3 = /3 C [1 + 5(U)}. 

III. OFF-EQUILIBRIUM SCALING THEORY AT 
FIRST-ORDER TRANSITIONS 

Analogously to the case of off-equilibrium slow dynam¬ 
ics at continuous transitions [ 2 , : 4|, we construct an off- 
equilibrium scaling theory in terms of the effective length- 
scale exponent v describing the finite-size scaling (FSS) 
at equilibrium, and the dynamic exponent 2 associated 
with the equilibrium relaxational dynamics correspond¬ 
ing to the heat-bath or Metropolis updating algorithm. 
In the following we mainly focus on the FOTs of 2D Potts 
models, but most scaling arguments can be straightfor¬ 
wardly extended to other FOTs. 


A. Equilibrium exponents 

1. The length-scale exponent v 


In the case of FOTs the effective length-scale exponent 
v controlling the FSS at Tf generally depends on the ge¬ 
ometry of the lattice |.38l4 l l| . i.e., whether it is square L 2 
or slab-like L\ x L\\ with L || Tj_, and on the boundary 
conditions fTot - PT^ . 

In the case of square L 2 systems, the FSS around T c is 
generally described by the effective length-scale exponent 
1 / = 1/d = 1/2 (38. IMt) . This length-scale exponent 
may be associated with a discontinuity fixed point in the 
renormalization-group framework [43]. It represents the 
limiting case of continuous transitions, leading to the en¬ 
ergy discontinuity at T c in the thermodynamic limit [ 44 J . 

In the case of slab-like geometries, the length-scale ex¬ 
ponent controlling the FSS at T c with res pec t to the fi¬ 
nite size L l may significantly change [33, [40|, [41J. This 
is essentially due to the anisotropic behavior of the lon¬ 
gitudinal length scale £||, which may show a nontrivial 
power-law (or exponential) dependence on Tj_, such as 
$11 ~ £j_- This may give rise to a change of the effective 
dimensions of the FSS in a slab at the FOT. Indeed, as¬ 
suming that the free-energy density scales as the inverse 
of the relevant critical volume V c ~ L± x £|| ~ 
we obtain the effective dimension d e g- = 1 + e. Thus, in 
this anisotropic setting, the effective length-scale expo¬ 
nent controlling the FSS with respect to L± turns out to 
be 


1 _ 1 
d e ff 1 + e 


(9) 


This value allows for the discontinuity of the energy den¬ 
sity at T c after taking the L± — > 00 limit. Eq. m can 
be also obtained from the FSS of the corresponding one- 
dimensional quantum model at the first-order quantum 
transition E(J [ij|, exploiting the quantum-to-classical 
mapping [48l ] which relates the gap A (energy difference 
of the lowest levels) of the quantum model to the in¬ 
verse longitudinal correlation length £|| of the 2D classical 
model in a slab geometry. 

In a slab geometry with symmetric boundary condi¬ 
tions, such as PBC, the relevant configurations involve 
domain walls which divide the system into successive re¬ 
gions of high-T and low-T phases, whose length scale £|| 
is expected to diverge exponentially with the transverse 
size Tj_. This length scale is related to the interfacial 
tension er, i.e., £|| ~ e aL . An analogous scenario applies 
to 2D Ising models in a slab geometry, at their FOTs 
driven by the magnetic field in their low-T phase [38}. 
As a consequence, Eq. © leads to the extreme value 
v —f 0. 

The situation changes in the case of MBC, for which 
increases as a power law of L±. In particular, in the case 
of the FOTs of 2D Potts models we have £|| ~ L, i.e., 
e = 1. This can be inferred from the behavior of the gap, 
A ~ 1/T, of the one-dimensional quantum Potts model 
at the transition with the corresponding MBC [4l| . 

Note that e = 1 is not a general feature of MBC. For 
example, at the low-T coexistence line of 2D Ising mod¬ 
els, where FOTs are driven by an external magnetic field 
coupled to the order parameter, we have £m ~ L 2 (thus 
e = 2) in the case of MBC favoring the two different mag¬ 
netized phases (i.e., fixed and opposite boundary condi¬ 
tions for the order-parameter field) [Ic|. 

In conclusion, 2D Potts models in square geometries 
and slab geometries with MBC share the same length- 
scale exponent 


v = 1/2. 


( 10 ) 


2. The equilibrium dynamic exponent z of the purely 
relaxational dynamics 

The equilibrium dynamic exponent z of the relax¬ 
ational dynamics is related to the equilibrium large-T 
behavior of the autocorrelation time r of observables at 
T c , i.e., t ~ L z . Again we expect that it depends on the 
boundary conditions at FOTs. 

In the case of symmetric boundary conditions, such 
as PBC, the autocorrelation time at T c is expected to 
exponentially increase with increasing lattice size, cor¬ 
responding to z 00 . This is related to the exponen¬ 
tial increase of the tunne ling time between the coexisting 
phases at FOTs, indeed (49| r ~ e aL (neglecting powers 
of L in the prefactor) where cr is the interfacial free en¬ 
ergy per unit length. To overcome this exponential slow 
down at FOTs, multicanonical updating algorithms have 
been developed j49| . 
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FIG. 1: (Color online) Log-log plot of the integrated auto¬ 
correlation time r of the magnetization at T c vs L, computed 
by equilibrium heat-bath MC simulations. The dotted and 
dashed lines show fits of the data for the lattice sizes L > 16, 
to the ansatzes r = oL 3 + bL 2 , corresponding to z = 3 with 
the expected 0(1/L) corrections, and to r = aL z . These 
fits are practically equivalent, indeed the corresponding lines 
are hardly distinguishable (both of them give an acceptable 
xVd.o.f. « 1). 



The behavior of the autocorrelation time drastically 
changes in the case of MBC. In this case the dynamics 
at T c is essentially related to the interface enforced by 
MBC, which moves along the slab. Indeed, this gives 
rise to a power-law behavior: r ~ L z . 

We numerically estimate z by equilibrium MC sim¬ 
ulations of the q = 20 Potts model at T c for slab-like 
geometries (in particular for anisotropic L± x L\\ with 
L± = 2L + 1 and L\\ L), with MBC. Some details on 

the calculation of the autocorrelation time of the magne¬ 
tization and energy are reported in App. [B] In Fig. |T] we 
show data for the integrated autocorrelation time r of the 
magnetization ([3]), from L = 12 to L = 32 and L || = 8 L 
(as already discussed above, L\\ ~ L is the correct scaling 
of the longitudinal size since £|| ~ L). 

These results strongly support a power-law behavior, 
i.e., t ~ L z . In order to estimate z , we fit the data to the 
ansatzes r = aL z [fitting tha data for L > 16, it gives z = 
3.08(2) and a = 0.33(2)] and r = aL 3 +bL 2 corresponding 
to z = 3 with the expected 0(1/L) corrections [it gives 
a ss 0.46(1) and b = —0.7(2)]. These fits work equally 
well as shown in Fig. [I] We consider 

2 = 3.0(1) (11) 

as our final estimate (the central estimate 2 = 3 is also 
supported by the off-equilibrium simulations, see below). 
The integrated autocorrelation time of the energy density 
gives substantially equivalent results. 

Although this estimate of z is obtained for heat-bath 
MC simulations at q = 20, we expect that it holds for 
any q > 4. Indeed, the value of z should be intrinsically 
related to the interface dynamics which is expected to be 
shared by all q > 4 at their FOTs. Moreover, it should 


extend to the whole class of purely relaxational dynam¬ 
ics (model A according to Ref. [33]), including also the 
Metropolis upgrading. 


B. Off-equilibrium scaling ansatzes 

A scaling theory for the off-equilibrium dynamics 
across T c can be heuristically derived by scaling argu¬ 
ments, similar to those commonly used at continuous 
transitions. 

Assuming the existence of a nontrivial scaling behavior 
around t = 0 corresponding to T c , we may expect that it 
is controlled by two scaling variables, such as 

n = (t/t s )~ v /L, r 2 = t/L z . (12) 

In particular n may be interpreted as the ratio between 
the equilibrium correlation length at the given /3(t) and 
L. The off-equilibrium scaling behavior arising from the 
protocol {8]) is meant to describe the deviations of the sta¬ 
tistical correlations from the equilibrium scaling behav¬ 
ior, when they cannot adapt themselves to the changes 
of the temperature across T c . Thus equilibrium scaling 
should be recovered in the limit |t "2 | —> oo keeping n 
fixed. 

Equivalently we consider the scaling variables 

u = tg/L, k = v/(l + zv), (13) 

w = K t = zk, (14) 

which are combinations of tt ad r 2 . 

We already note that the case of symmetric boundary 
conditions, such as PBC, appears problematic, due to 
the divergence of the dynamic exponent z. In particular, 
taking the z —> oo limit of the exponents k and n t in the 
case of a square geometry, one would naively obtain k = 0 
and Kt = 1. This may hint at the absence of a nontrivial 
scaling behavior around T c , i.e., we may not observe a 
nontrivial off-equilibrium scaling behavior around t = 0 
in the off-equilibrium protocol (]5]|. Alternatively, this 
may indicate a logarithmic scaling behavior, for example 
with scaling variables u ~ In (t s )/L and w ~ t/t s . As we 
shall see, numerical simulations favor a regular behav¬ 
ior around t = 0 extending to t/t s > 0, which may be 
somehow related to a metastability phenomenon. 

Systems with MBC appear more promising to realize 
an off-equilibrium scaling behavior. Indeed the corre¬ 
sponding values of the equilibrium exponent v and z pro¬ 
vide well defined exponents k and Kt in Eqs. m and 
(flTl) . In the case of the FOT of the 2D Potts model, for 
which v = 1/2 and z = 3.0(1), we obtain k = 0.200(4) 
and k t = 0.600(8), which lead to a nontrivial power-law 
dependence of the scaling variables. Therefore, observing 
an off-equilibrium scaling behavior around t = 0 is to be 
expected in this case. 

Our main working hypothesis is that the slow dynamics 
across T c presents a double scaling behavior in the large- 
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L limit, in terms of the scaling variables u and w, i.e., 

m r (t,t s ,L) S3 f m (u,w), (15) 

e r (t, t s , L) S3 f e {u, w). (16) 

We also expect that, if we start from the high-T phase 
[sign + in Eq. ( 0 ], the following asymptotic limits apply: 

lirn f#(u,w)= 0, lim f#(u, w) = 1, (17) 

w—>—oo w—t oo 

corresponding to the large-T equilibrium values at the 
two phases (the subscript # corresponds to both m and 
e). The limits are reversed if we start from the low-T 
phase [sign - in Eq. 0]. 

Summarizing, in the cases where an off-equilibrium 
scaling behavior is driven by slow variations across T c , 
around t = 0 corresponding to T(t) = T c , Eqs. (fl5l) 
and (Hsu with the scaling variables m and m are ex¬ 
pected to provide the asymptotic scaling behaviors, when 
L is much larger than any other length scale, thus when 
L ^ (see the end of Sec. Ill AD . These asymptotic 
behaviors should be approached with power-law 0 (L _ “) 
suppressed corrections, presumably u> = 1. Note that 
for observables depending on the spatial coordinates one 
may add the spatial scaling dependence on x/L. 

In the case of MBC the off-equilibrium behavior is ex¬ 
pected to be related to the dynamics of the interface en¬ 
forced by the MBC. This scenario implies a close relation 
between the scaling functions f m and f e of the magneti¬ 
zation and energy density. Let us assume that the slowest 
modes, determining the off-equilibrium dynamics of the 
magnetization and energy density, are associated with 
the wall separating the spatial regions corresponding to 
the two different phases. If we define x(t) as the loca¬ 
tion of the wall along the x\ axis, separating the low-T 
from the high-T region, we expect that m r = e r = 1 for 
x\ < x(t) and m r = e r = 0 for X\ > x(t). Thus, when 
the width of the domain wall becomes negligible, thus 
asymptotically for large L , we expect that 

m r (t) ss e r (t) ss ^[1 + X w (i)], (18) 

where X w (t) = x(t)/L with —1 < X w (t) < 1. As a trivial 
consequence, we would have 

fe(u,w) = fm(u,w). (19) 

Therefore, X w (u,w) = 2 f m (u,w) — 1 may be considered 
as an estimator of the average position of the interface in 
the large-T limit. 

The off-equilibrium scaling behaviors (I13II16I) at FOTs 
are analogous to those expected at continuous transitions 
when slowly crossing the critical temperature, usually re¬ 
lated to the KZ mechanism Hi- The main difference 
is related to the fact that at continuous transitions the 
critical exponents, such as v and z, do not depend on the 
boundary conditions (only the scaling functions do). For 
example, in the case of the 2D Ising model corresponding 


to the q = 2 Potts model, the off-equilibrium dynamics 
of the magnetization 0 across the continuous transition 
is expected to be 

m{t, t s , L ) = L~ v/2 f m (u , w) (20) 

where 77 = 1/4 and the scaling variables have the same 
form as those in Eqs. m and ED- The exponents k 
and Kt are obtained using the Ising exponents v = 1 
and z ss 2.167 associated with the purely relaxational 
dynamics (see, e.g., Refs. [50,.51] and references therein), 
i.e., k ss 0.316 and Kt ~ 0.684. At the Ising critical point, 
different geometries and boundary conditions can only 
change the scaling function f m . 

We finally mention that similar scaling arguments have 
been also employed to study the effects of smooth spa¬ 
tial inhomogeneities at FOTs [52l - l54l |. for example, when 
we assume a spatially dependent temperature and look 
at the behavior of the system around the spatial region 
where the temperature takes the value at the transition 
point of the homogenous system. 

IV. MONTE CARLO SIMULATIONS 

In order to check the predictions of the previous sec¬ 
tion, we present a numerical analysis of MC simulations 
following the off-equilibrium protocol 0 , for various ge¬ 
ometries, boundary conditions, lattice sizes and time 
scales. We use the heat-bath upgrading, see App. [A] 
varying the temperature according to Eq. 0 , after each 
sweep which corresponds to a unit time. 

A. Square systems with PBC 

To begin with, we present results for the square L 2 lat¬ 
tice with PBC, and the off-equilibrium protocol 0 start¬ 
ing from the high-T phase. In this case the magnetization 
vanishes by symmetry, as a consequence of the average 
on the initial Gibbs ensemble at /3(t*) < 0. Therefore 
we look at the behavior of the energy density. Figure [2] 
shows some results for the renormalized energy density, 
cf. Eq. 0 , obtained by heat-bath MC simulations of the 
q = 20 Potts model for various L and time scales t s . 

When increasing t s and T, the data approach the cor¬ 
responding L —> 00 equilibrium values up to t = 0, 
where e r vanishes within errors (we recall that e r = 0 
is the T —> T+ limit of the equilibrium value in the 
thermodynamic limit). Then they increase slowly up 
to t/t s = t* > 0 with r* ss 0.005, corresponding to 
T(t) S 3 0.995 T c , remaining well below the equilibrium 
values of the low-T phase. Then the data show a sharp 
crossover to the values corresponding to the low-T phase. 
Note also that the lower panel of Fig. [2] reports data with 
In (t s )/L si const, which do not support scaling with re¬ 
spect to the scaling variables u ~ hi(t s )/L and w ~ t/t s . 
Of course, we cannot exclude that an eventual logarith¬ 
mic scaling behavior may set in for larger L and t s . 
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FIG. 2: (Color online) Data of e r {t) for the 2D q = 20 Potts 
model in square L 2 lattices with PBC, starting from the high- 
T phase. We show two sets of data: (top) data for large 
lattices with t s = O(10 6 ); (bottom) data for smaller lattices, 
for which the ratio ln(t s /16 )/L « 0.2558 is kept fixed (to 
check the possibility of a logarithmic scaling around t = 0). 
Both of them indicate a smooth behavior up to t/t s « 0.005 
and then a sharp relaxation toward the low-T values, without 
showing nontrivial scaling behaviors around t = 0. 

The numerical results for medium-size L = 0{ 10 2 ) lat¬ 
tices and time scales t s = O(10 6 ) suggest a smooth non¬ 
singular behavior around t/t s = 0, without hinting at 
nontrivial off-equilibrium scaling behaviors around t = 0. 
This scenario was somehow anticipated in the previous 
section, as a consequence of the trivial values of the ex¬ 
ponents k and Kt obtained using the scaling ansatzes m 
and (HU). It may be related to some form of metastabil¬ 
ity developing in this slow cooling procedure, which likely 
requires another theoretical framework. This issue calls 
for further investigation. 

We expect that analogous scenarios occur at the FOTs 
of any q > 4. 


B. Results for slab-like systems with MBC 

We now present results for the 2D Potts models with 
q = 10 and q = 20, in anisotropic (2 L + 1) x L\\ lattices 


FIG. 3: (Color online) Data of e r {t) (top) and m r (t ) (bottom) 
for the 2D q = 20 Potts model with MBC, in anisotropic 
L± x L || lattices with L± = 2L + 1 and L« = 8 L, for L = 16 
and L = 32, and various values of t 3 . 
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FIG. 4: (Color online) Data of m r (t = 0) for q = 20, start¬ 
ing from high-T (he) and low-T (ch) phases, and for various 
values of u = t° s 2 /L (which is kept fixed when varying t 3 and 
L ) vs 1/L. They appear to approach a constant value with 
increasing L, supporting the general scaling ansatz 03 with 
k = 0.2. The dotted lines show a linear fit of the data for the 
largest available lattice sizes. 

























7 


with MBC, starting from the high-T and low-T phase. 
We consider the slab-limit L y Lj_ , for which numerical 
results can be straightforwardly obtained by increasing 
the longitudinal size up to the point where the data ap¬ 
pear stable within the errors. We checked that Ty = 8 L 
turns out to be sufficiently large to effectively provide 
infinite-Tii results within the errors (the linear scaling of 
L || with L takes into account that ~ L± for MBC). In 
our numerical simulations we choose |<$(tj)| = |<5(t/)| = 
1/32. However, as we shall see, the emerging scaling 
behavior does not depend on these particular values, be¬ 
cause it is essentially related to the behavior around t = 0 
where T(t) = T c . 

In Fig. [3] we report some raw data during a slow vari¬ 
ation of the temperature across T c for the slab geometry 
with MBC. They are obtained starting from the high-T 
phase and show the behavior of the energy and the mag¬ 
netization when the transition is crossed with different 
time scales t s and lattice sizes L. The data show that 
the effective passage from one phase to the other occurs 
around t/t s = 0. As we shall see, the dependence of 
the various curves on L and t s can be cast in the off- 
equilibrium scaling behavior put forward in the previous 
section, cf. Eqs. (11311161) . 

In order to check the scaling in the variable u = t*/L, 
we first note that Eqs. (usd and m imply that at t = 0 

m r (0,t s ,L) ss g m (u), e r (0, t s , L) ss g e (u). (21) 

Therefore, we expect that data at t = 0 and fixed 
u = t*/L must converge to nontrivial w-dependent values 
with increasing L. Fig. [I] shows data at some fixed val¬ 
ues of u (using k = 0.2 obtained taking the central value 
2 = 3), for q = 20 and q = 10. They appear to converge 
to nontrivial values, supporting the above asymptotic be¬ 
havior, with corrections which approximately decay as 
OiL- 1 ). 

Analogous results are obtained by slightly changing the 
value of z, according to the equilibrium estimate 2 = 
3.0(1), corresponding to k = 0.200(4). Actually, one may 
assume the off-equilibrium scaling m to independently 
estimate 2 from the off-equilibrium data. For example, 
by allowing for deviations from 2 = 3, i.e., 2 = 3 + ^ 2 , 
and interpreting the scaling corrections of the data of 
Fig. |4] as due to Sz, we find again that the optimal value 
is 2 = 3 with a few percent of uncertainty (corrections to 
scaling have different sign in some cases, which cannot 
be explained by a unique shift of 2 , thus 2 ~ 3 appears 
as the optimal value). 

Note also that, since the equilibrium average position 
of the interface at T c is expected at equal distances from 
the boundaries in the large- L limit, leading to the asymp¬ 
totic equilibrium values m r (T c ) = e r (T c ) = 1/2, we ex¬ 
pect that lim „_ > , 00 g#{u) = 1 / 2 . 

The scaling with respect to w at fixed values of u 
is supported by the plots in Figs. [3 [ 6 j and [7] respec¬ 
tively for q = 20 with hot and cold starting point and 
q = 10 with hot start. In all cases the data approach 
an asymptotic function of the scaling variable w, as pre- 



FIG. 5: (Color online) Data of m r (t) and e r {t) for q = 20 
starting from the high-T phase, i.e., T > T c , for two values of 
u = tg/L with re = 0.2, i.e., u ss 0.7578 (top) and u « 0.5743 
(bottom). With increasing L, they approach asymptotic 
curves when plotted versus w = t/tg* with ret = 0.6, sup¬ 
porting the scaling ansatzes m and m ■ Moreover, the 
data are consistent with the relation G3 predicting the same 
asymptotic curve for m r and e r . 



FIG. 6: (Color online) Data for m r (t) and e r (t) starting from 
the cold phase, i.e., T < T c , keeping the scaling variable u tn 
0.5743 fixed. They approach a unique asymptotic curve with 
increasing L, when plotted versus w = t/t £*, supporting the 
scaling ansatzes m and m, and the interface relation m ■ 











FIG. 7: (Color online) Data of m r (t) for q = 10 versus w = 
t/t? at a fixed value of u = tl/L as 0.3789 (with Kt = 0.6 
and re = 0.2), starting from the high-T phase. They approach 
an asymptotic scaling curve, in agreement with Eq. (USD- 



FIG. 8: (Color online) Check of universality between the 
heat-bath and Metropolis dynamics. The scaling curves coin¬ 
cide after a rescaling of the scaling variables u and w of the 
Metropolis data. In particular, the rescaling of the variable u 
can be fixed by looking for data with equal m r { 0) (we obtain 
that u as 0.5743 of the heat-bath dynamics approximately 
corresponds to u = 1 of the Metropolis dynamics). 

dieted by the scaling theory of Sec. IIIII Moreover, the 
data of m r (t) and e r (t) shown in Figs. [5] and [6] approach 
the same scaling curve, in agreement with the asymptotic 
relation (USD, thus fully supporting the hypothesis that 
the off-equilibrium behavior is essentially controlled by 
the time-dependent position of the interface. 

Finally, in Fig. [8] we compare results obtained using 
the heat-bath algorithm with those obtained by a stan¬ 
dard Metropolis algorithm with only one trial per site, 
see App. [A] The curves match after a trivial rescaling of 
the scaling variables u and w, supporting the expected 
universality with respect to the type of relaxational dy¬ 
namics. 

We conclude that our numerical analysis for slab ge¬ 
ometries with MBC supports the off-equilibrium scaling 
behaviors of the magnetization and energy density put 


forward in Sec. iotbi cf. Eqs mm- 


V. CONCLUSIONS 

We investigate off-equilibrium behaviors at FOTs 
driven by a time dependence of the temperature across 
the transition point T c . Usually, off-equilibrium behav¬ 
iors at FOTs are associated with phenomena of metasta¬ 
bility and hysteresis [46] . We focus on the possibility 
of nontrivial off-equilibrium scaling behaviors driven by 
slow changes of the temperature, analogous to those aris¬ 
ing at continuous transitions, leading to the KZ mech¬ 
anism 0 i- When slowly varying the temperature 
across a continuous transition, for example linearly as 
T(t)/T c w 1 - t/t s , the KZ mechanism predicts a non¬ 
trivial off-equilibrium scaling behavior in the limit of 
slow variations, due to the diverging length scale at 
T c . We show that phenomena analogous to the KZ off- 
equilibrium scaling emerge also at FOTs, when appropri¬ 
ate boundary conditions are considered. 

We consider the 2D Potts models, which provide an 
ideal testing ground to investigate issues related to FOTs. 
In our discussion we consider a purely relaxational dy¬ 
namics such as that obtained by heat-bath and Metropo¬ 
lis upgrading in MC simulations. We study the off- 
equilibrium behavior in the case of a time-dependent 
temperature crossing the FOT. In particular, we con¬ 
sider a linear dependence of the inverse temperature 
/3(f) = /3 C (1 ± t/t s ), starting from the high-T or low-T 
phase. 

We point out that off-equilibrium behaviors at FOTs 
are extremely sensitive to the geometry and boundary 
conditions of the system. This peculiar dependence is es¬ 
sentially related to the equilibrium relaxational dynamics 
at T c . For symmetric boundary conditions, such as PBC, 
we expect an exponentially slow dynamics due to an ex¬ 
ponentially large tunneling time r ~ e aL . On the other 
hand, a power-law slowing down r ~ L z with z ss 3.0 
is found when considering MBC, i.e., when the bound¬ 
ary conditions at two opposite sides of the system are 
related to the different high-T and low-T phases, effec¬ 
tively generating an interface separating the coexisting 
phases. We argue that an off-equilibrium scaling behav¬ 
ior around t = 0 (i.e., the time corresponding to T c ) is 
realized for MBC. This is controlled by the correspond¬ 
ing equilibrium length-scale and dynamic exponents. We 
argue that this scaling behavior is essentially related to 
the dynamics of the interface enforced by MBC. 

In the case of slab-like geometries with MBC, the nu¬ 
merical results for the q = 10 and q = 20 Potts model 
support the emergence of an off-equilibrium scaling pic¬ 
ture characterized by power-law behaviors, analogous to 
those predicted by the KZ theory at continuous transi¬ 
tions. 

On the other hand, symmetric boundary conditions 
do not apparently lead to nonanalytic scaling behaviors 
at t = 0, but rather to a delayed sharp relaxation to 
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the other phase at t/t s > r* with r* > 0, which may 
be related to a metastability phenomenon arising from 
the slow dynamics across the FOT. This point deserves 
further investigation to physically understand it. Fur¬ 
ther checks of the observed asymptotic behavior may also 
be called for. Indeed, we cannot exclude the possibility 
that a different (logarithmic) asymptotic behavior sets 
in for sizes and time scales larger than those considered 
in our numerical analysis, which are L = O(10 2 ) and 

t, = o(io 6 ). 

Our scaling arguments are quite general, therefore they 
should also apply to higher-dimensional systems, such as 
the 3D Potts models that undergo FOTs. Such an exten¬ 
sion may generally depend on the geometry of the sys¬ 
tem, e.g., cubic-like L 3 , slab-like Lj_ x L jj -1 and tube-like 

L 2 ± x L|| (with L || L± ) geometries, as well as on the 
boundary conditions. In particular, we again expect that 
geometries and boundary conditions favoring the emer¬ 
gence of an interface should give rise to off-critical scaling 
behaviors similar to those of the KZ theory at continuous 
transitions. This issue calls for further investigation. 

Off-equilibrium scaling behaviors may also appear at 
FOTs driven by magnetic fields. For example, one may 
consider 0(fV)-symmetric spin models in the low-T or¬ 
dered phase, where FOTs are driven by an external mag¬ 
netic field coupled to the spin variables. Then, one may 
consider the off-equilibrium dynamics driven by a time- 
dependent magnetic field h(t) = t/t s across the transition 
point h = 0. We expect that in the case of Ising models 
(N = 1) an off-equilibrium scaling behavior may emerge 
in systems with boundary conditions enforcing the pres¬ 
ence of an interface, analogously to what is observed at 
the thermal FOTs of the 2D Potts models. In the case 
of a continuous symmetry (N >2), off-equilibrium scal¬ 
ing behavior may emerge from the spin-wave dynamics 
(Goldstone modes related to the broken O (N) symme¬ 
try) [39]. We believe that these issues are worth being 
further investigated. 

One may also consider evolutions different from the 
purely relaxational ones. They would generally lead to 
other values of the dynamic exponent z. 

We also mention that off-equilibrium behaviors arising 
from sudden quenches below and at the transition point 
have been discussed in several works, see, e.g., [55l - l6l| . 

Off-equilibrium behaviors at FOTs are quite general, 
they should be observable in many physical contexts 
where the FOTs are approached by varying the system 
parameters. The off-equilibrium protocol investigated in 
this paper may be exploited to probe the main features of 
systems at the FOT. Moreover, our results may turn out 
useful in understanding more complicated off-equilibrium 
phenomena at FOTs. For example, as a case of physical 
interest we mention the effects of the intrinsic space-time 
inhonrogeneity of the quark-gluon plasma formation in 
heavy-ion collisions [62], whose equilibrium T-/z (/x is the 
chemical potential) phase diagram is expected to have 
a FOT line [63 which may be crossed during heavy-ion 
collisions. Another interesting context is that of the uni¬ 


verse cosmology, jwhich was the original ground of the 
Kibble proposal p| to understand the effects of an ex¬ 
panding universe through a continuous transition. From 
analogous studies of off-equilibrium behaviors at FOTs 
we may learn the effects of an expanding and cooling 
universe that passes through a FOT. 

Analogous off-equilibrium phenomena should be also 
observable in quantum many-body systems, at first- 
order quantum transitions. Some issues arising from 
slow (quasi-adiabatic) passages through quantum FOTs 
have been recently discussed, in particular for some one- 
dimensional quantum chains [42], [HI, [HI, including is¬ 
sues related to adiabatic evolutions in quantum compu¬ 
tations [il. l66i - [68j . 
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Appendix A: Metropolis and Heat-bath updatings of 
the Potts models 

A heat-bath updating of a single site variable consists 
in the change s x —> s' x with probability ~ e ~ H ( s x )/ T 
independent of the original spin s x . 

The Metropolis updating of a single spin s x is per¬ 
formed by (i) proposing a new spin ^ s x by taking one 
of the other q— 1 states with equal probability, (ii) accept¬ 
ing the change with probability Min[e[' ff ^ Sx ^ _i? ^ Sx ^/ r , 1]. 

The heat-bath updating is generally more effective 
than a single Metropolis updating, because the new 
variable is not correlated with the previous one. The 
Metropolis updating tends to be equivalent to the heat- 
bath one when a large number of trials are performed. 

Both updating procedures give rise to a purely relax¬ 
ational dynamics, usually named model A [37[ . whose 
class also includes configuration updatings by Langevin 
equations with white noise. The time unit during the re¬ 
laxational dynamics is generally associated with a com¬ 
plete sweep of the lattice variables. 


Appendix B: Computation of the equilibrium 
autocorrelation time 

The integrated autocorrelation time of a given quantity 
Q is defined as 


_ 1 C(t) 

T ~ 2 ^ C( 0)’ 

t= — OO v ' 


(Bl) 


where C(t ) = (( Q(t ) — (Q)) (Q(0) — (Q))) is the autocor¬ 
relation function of Q (t is the discrete Monte Carlo time, 
where a time unit is given by a sweep, i.e., a heat-bath 
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update of all lattice variables). Averages are taken at 
equilibrium. Estimates of the corresponding integrated 
autocorrelation time r can be obtained by the binning 
method (see, e.g., Ref. [el Iz3 for discussions of this 
method and its systematic errors), using the estimator 


where Eq is the naive error calculated without taking into 
account the autocorrelations, and E is the correct error 
found after binning, i.e., when the error estimate becomes 
stable with respect to increasing of the block size b. The 


statistical error At is just given by Ar/r = \J 2 /rib where 
rib is the number of blocks corresponding to the estimate 
of E. As discussed in Ref. [69[ this procedure leads to 
a systematic error of 0{r/b). In our cases the ratio r/b 
will always be much smaller than the statistical error, so 
we will neglect it. Eq. (1B2I) can be easily extended to the 
case the quantity Q is measured every n m sweeps, i.e., 
t = n m E 2 /(2 Eq), which is of course meaningful only if 
n m < r. 

Of course, r depends on the quantity Q considered. 
However, barring unlikely exceptions, all quantities lead 
to the same asymptotic power-law behavior r ~ L z . 


[1] T.W.B. Kibble, J. Phys. A 9, 1387 (1976). 

[2] W.H. Zurek, Nature 317, 505 (1985). 

[3] S. Gong, F. Zhong, X. Huang, and S. Fan, New J. Phys. 
12, 043036 (2010). 

[4] A. Chandran, A. Erez, S.S. Gubser, and S.L. Sondhi, 
Phys. Rev. B 86, 064304 (2012). 

[5] W.H. Zurek, U. Dorner, and P. Zoller, Phys. Rev. Lett. 
95, 105701 (2005); J. Dziarmaga, Phys. Rev. Lett. 
95, 245701 (2005); A. Polkovnikov, Phys. Rev. B 72, 
161201(R) (2005). 

[6] A. Polkovnikov and V. Gritsev, Nature Phys. 4, 477 
(2008). 

[7] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalat- 
tore, Rev. Mod. Phys. 83, 863 (2011). 

[8] I. Chuang, R. Durrer, N. Turok, and B. Yurke, Science 
251, 1336 (1991). 

[9] M. J. Bowick, L. Chandar, E. A. Schiff, A. M. Srivastava, 
Science 263, 943 (1994). 

[10] C. Bauerle, Yu M. Bunkov, S. N. Fisher, H. Godfrin, and 
G. R. Pickett, Nature 382, 332 (1996). 

[11] V. M. H. Ruutu, V. B. Eltsov, A. J. Gill, T. W. B. Kibble, 
M. Krusius, Y. G. Makhlin, B. Placais, G. E. Volovik, 
and W. Xu, Nature 382, 334 (1996). 

[12] R. Carmi, E. Polturak, G. Koren, Phys. Rev. Lett. 84, 
4966 (2000). 

[13] S. Casado, W. Gonzalez-Vinas, H.L. Mancini, and S. 
Boccaletti, Phys. Rev. E 63, 57301 (2001). 

[14] R. Monaco, J. Mygind, and R. Rivers, Phys. Rev. Lett. 
89, 080603 (2002). 

[15] A. Maniv, E. Polturak, G. Koren, Phys. Rev. Lett. 91, 
197001 (2003). 

[16] S. Casado, W. Gonzales-Vinas, and H. Mancini, Phys. 
Rev. E 74, 047101 (2006). 

[17] R. Monaco, J. Mygind, M. Aaroe, J. Rivers, and V.P. 
Koshelets, Phys. Rev. Lett. 96, 180604 (2006). 

[18] L.E. Sadler, J.M.Higbie, S.R. Leslie, M. Vengalattore, 
and D.M. Stamper-Kurn, Nature 443, 312 (2006). 

[19] C.N. Weiler, T. W. Neely, D. R. Scherer, A. S. Bradley, 
M. J. Davis and B. P. Anderson, Nature 455, 948 (2008). 

[20] D. Golubchik, E. Polturak, and G. Koren, Phys. Rev. 
Lett. 104, 247002 (2010). 

[21] D. Chen, M. White, C. Borries, and B. DeMarco, Phys. 
Rev. Lett. 106, 235304 (2011). 

[22] S. C. Chae, N. Lee, Y. Horibe, M. Tanimura, S. Mori, B. 
Gao, S. Carr, and S.-W. Cheong Phys. Rev. Lett. 108, 


167603 (2012). 

[23] M.A. Miranda, J. Burguete, H. Mancini, and W. 
Gonzalez-Vinas, Pliys. Rev. E 87, 032902 (2013). 

[24] S. Eitemaee and P. C. Haljan, Phys. Rev. A 87, 
051401(R) (2013). 

[25] S. Ulm, S. J. Ronagel, G. Jacob, C. Degrinther, S. T. 
Dawkins, U. G. Poschinger, R. Nigmatullin, A. Retzker, 
M. B. Plenio, F. Schmidt-Kaler, K. Singer, Nat. Com- 
mun. 4, 2290 (2013). 

[26] K. Pyka , J. Keller, H. L. Partner, R. Nigmatullin, T. 
Burgermeister, D. M. Meier, K. Kuhlmann, A. Retzker, 
M. B. Plenio, W. H. Zurek, A. del Campo, and T. E. 
Mehlsaubler, Nat. Commun. 4, 2291 (2013). 

[27] G. Lamporesi, S. Donadello, S. Serafini, F. Dalfovo, and 
G. Ferrari, Nat. Phys. 9, 656 (2013). 

[28] L. Corman, L. Chomaz, T. Bienaime, R. Desbuquois, C. 
Weitenberg, S. Nascimbne, J. Dalibard, and J. Beugnon, 
Phys. Rev. Lett. 113, 135302 (2014). 

[29] N. Navon, A. L. Gaunt, R. P. Smith, and Z. Hadzibabic, 
Science 347, 167 (2015). 

[30] S. Braun, M. Friesdorf, S.S. Hodgman, M. Schreiber, J.P. 
Ronzheimer, A. Riera, M. del Rey, I. Bloch, J. Eisert, and 
U. Schneider, PNAS 112, 3641 (2015). 

[31] R.J. Baxter, Exactly solved models in statistical mechan¬ 
ics , (Academic Press, 1982). 

[32] F.Y. Wu, Rev. Mod. Phys. 64, 235 (1982). 

[33] A. Kliimper, A. Schadschneider, and J. Zittartz, Z. Phys. 
B 76, 247 (1989). 

[34] E. Buddenoir and S. Wallon, J. Phys. A 26, 3045 (1993). 

[35] F. Igloi and E. Carlon, Phys. Rev. B 59, 3783 (1999). 

[36] W. Janke and S. Kappler, Europhys. Lett. 31, 345 
(1995). 

[37] P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49, 
435 (1977). 

[38] V. Privman and M. E. Fisher, J. Stat. Phys. 33, 385 
(1983). 

[39] M. Fisher and V. Privman, Phys. Rev. B 32, 447 (1985). 

[40] M. Campostrini, J. Nespolo, A. Pelissetto, and E. Vicari, 
Phys. Rev. Lett. 113, 070402 (2014). 

[41] M. Campostrini, J. Nespolo, A. Pelissetto, and E. Vicari, 
Phys. Rev. E 91, 052103 (2015). 

[42] M. Campostrini, A. Pelissetto, and E. Vicari, Phys. Rev. 
E 91, 042123 (2015); arXiv: 1507.08391 

[43] B. Nienhuis and M. Nauenberg, Phys. Rev. Lett. 35, 477 
(1975). 




11 


[44] M.E. Fisher and A.N. Berker, Phys. Rev. B 26, 2507 
(1982). 

[45] M.S.S. Challa, D.P. Landau, and K. Binder, Phys. Rev. 
B 34, 1841 (1986). 

[46] K. Binder, Rep. Prog. Phys. 50, 783 (1987). 

[47] K. Vollmayr, J.D. Reger, M. Scheucher, and K. Binder, 
Z. Phys. B 91, 113 (1993). 

[48] S. Sachdev, Quantum Phase Transitions, (Cambridge 
Univ. Press, 1999). 

[49] B. A. Berg and T. Neuhaus, Phys. Rev. Lett. 68, 9 
(1992). 

[50] M. P. Nightingale and H.W.J. Blote, Phys. Rev. B 62, 
1089 (2000). 

]51] G. Costagliola and E. Vicari, J. Stat. Mech.: Theory Exp. 
(2011) L08001. 

]52] C. Bonati, M. D’Elia, and E. Vicari, Phys. Rev. E 89, 
062132 (2014). 

[53] M. Campostrini, J. Nespolo, A. Pelissetto, and E. Vicari, 
Phys. Rev. E 91, 022108 (2015). 

[54] H. Weissmann, N. M. Shnerb, and D. A. Kessler, 
arXiv: 1508.03453 

[55] D. Gunton, M. San Miguel, and P.S. Sahni, Phase Tran¬ 
sitions and Critical Phenomena, edited by C. Domb and 
J.L. Lebowitz, Vol. 8, p. 267 (Academic Press, London, 
1983). 

[56] A.J. Bray, Adv. Phys. 43, 357 (1994). 


[57] L. Schiilke and B. Zheng, Phys. Rev. E 62, 7482 (2000). 

[58] Y. Ozeki, K. Kasono, N. Ito, and S. Miyashita, Physica 
A 321, 271 (2003). 

[59] L. Granasy, T. Pusztai, and J.A. Warren, J. Phys.: Con- 
dens. Matter 16, R1205 (2004). 

[60] E. Lorenz and W. Janke, Europhys. Lett. 77, 10003 
(2007). 

[61] M. Pleimling and F. Igloi, Europhys. Lett. 79, 56002 
(2007). 

[62] J.D. Bjorken, Phys. Rev. D 27, 140 (1983). 

[63] K. Rajagopal and F. Wilczek, in At the Frontier of Parti¬ 
cle Physics / Handbook of QCD, M. Shifman, ed., (World 
Scientific); arXiv:hep-ph/0011333 

[64] S. Sukuki and A. Dutta, arXiv: 1503.01512 

65 A.P. Young, S. Knysh, and V.N. Smelyanskiy, Phys. Rev. 
Lett. 104, 020502 (2010). 

[66] M.H.S. Amin and V. Choi, Phys. Rev. A 80, 062326 
(2009). 

[67] S. Miyashita, H. Do Raedt, and B. Barbara, Phys. Rev. 
B 79, 104422 (2009). 

[68] C. R. Laumann, R. Moessner, A. Scardicchio, and S. L. 
Sondhi, Phys. Rev. Lett. 109, 030502 (2012). 

[69] U. Wolff, Comput. Phys. Commun. 156, 143 (2004). 

[70] L. Del Debbio, G. M. Manca, and E. Vicari, Phys. Lett. 
B 594, 315 (2004). 


